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Abstract - We probe the diffusive motion of particles in slowly sheared three dimensional gran- 
ular suspensions. For sufficiently large strains, the particle dynamics exhibits diffusive Gaussian 
statistics, with the diffusivity proportional to the local strain rate — consistent with a local, quasi 
static picture. Surprisingly, the diffusivity is also inversely proportional to the depth of the parti- 
cles within the flow — at the free surface, diffusivity is thus ill defined. We find that the crossover 
to Gaussian displacement statistics is governed by the same depth dependence, evidencing a non- 
trivial strain scale in three dimensional granular flows. 



Flowing soft disordered solids, such as colloidal glasses 
and g els (l} |3], foams and emulsions [4jj6] or granular me- 
dia |7|- |13| , exhibit rich particle dynamics and complex rhe- 
ology. In the last decade, major progress has been made in 
the understanding of the rheology of granular media and 
suspensions. Important developments include the descrip- 
tion of rapid, "inertia!" dry granular flows [TJ[8] , and the 
observation of connections between this inertial rheology 
and the classical rheology of density matched suspensions 



14 . The most challenging regime is that of very slow 



flows, where the stresses become (approximately) rate in- 



dependent 15 and nonlocal effects play an important role 



[16 ■ 18 . Little is known about the fluctuations of the mi- 
croscopic grain motion, i.e., self-diffusion, in such flows. A 
recurring problem is that granular media are opaque, so 
only motion in two dimensional (2D) model systems [jjj, 
near transparant walls 10 19 or at a free surface can be 
observed 11 . Walls lead to layering [20] and other arti- 



facts 21 , while at the free surface, the particles experience 



a different local geometry than in the bulk. Measurements 
in the bulk of a granular flow are thus essential. 



Recently, we have shown 13 that for sufficiently slow 



flows, dry granular media and non-density, but optically 
matched suspensions exhibit identical flow patterns and 
rheology, opening up the possibility to probe 3D micro- 
scopic particle dynamics within slow granular flows. Here 



we experimentally probe the full 3D particle trajectories 
of a sheared granular suspension in a split bottom cell 
(Fig. 1). We focuss on sufficiently slow strain rates so 
that viscous effects are negligible: as detailed below, typ- 
ical Reynolds and Stokes numbers are of order 10 -4 and 
10, respectively, implying that we are in the so-called free 



fall regime, relevant for dry granular flows 13 22 



In a quasistatic picture, the mean squared displace- 
ments should only depend on the strain, or equivalcntly, 
the diffusion coefficient should be proportional to the 
strain rate, as seen already in [9j. By extracting the par- 
ticle diffusivity and local strain rates from our particle 
trajectories, we confirm that this is true in 3D also. 

The crucial point is that our 3D geometry allows us 
to probe the depth dependence of the diffusivity. In the 
simplest picture, one might expect the diffusion to the in- 
dependent of depth. However, our main result is showing 
that the diffusion is inverse proportional to the depth, and 
thus diverges near the free surface. To answer whether this 
depth dependence stems from a depth dependence of the 
mean squared displacements or of the characteristic strain, 
both of which can lead to a depth dependent diffusivity, we 
probe the displacement distributions. These evolve from 
very wide distributions at short time with (reduced) kur- 
tosis exceeding 40, to Gaussian distributions at late times. 
We find that, for a given depth, the kurtosis scales as a 
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Fig. 1: (Color online) (a) Setup: a slowly rotating disk at 
the bottom of a cubic container drives a slow granular flow. In 
this experiment, we probe driving rates Q, from 0.005 to 0.05 
rpm. The container is filled with non-density matched sus- 
pension consisting of PMMA beads, index matched fluid and 
fluorescent dye. The motion of the particles in this granular 
suspension are captured in 3D by vertical scanning of a laser- 
sheet, (b) Example of image obtained by illuminating a single 
slice and capturing the image with a CCD camera, (c) After 
particle tracking, the motion of the particles can be resolved in 
3D — in this rendered image, the color of the particles repre- 
sents their instantaneous velocity, (d) Strain rate field (for the 
case Q = 0.05 rpm), where the color code denotes the strain 
rate j. Symbols indicate the region in which we will probe local 
strain rates and fluctuations (circles: data used in Fig. 2abc; 
squares: additional data used in Fig. 2d). 



powerlaw with strain, and that we can collapse kurtosis 
data taken at different depths by introducing a depth de- 
pendent characteristic strain scale — the same depth de- 
pendent characteristic strain scale governs the diffusivity. 
Taken together, the granular fluctuations provide strong 
evidence for a non-trivial, depth dependent strain scale 
which governs slow granular flows. 

Setup — We probe the particle trajectories in the bulk of 
a gravitational (non-density matched) granular suspension 

This 



sheared in the split-bottom geometry 12 13 23 



geometry produces smooth and reproducible flow profiles, 
is well studied and allows to investigate a wide range of 
shear rates within a single experiment. This allows us to 
probe the relation between particle displacements, 7, and 
distance to the free surface. 

We use an identical geometry as the one used in 13 
(see Fig. la). A cubic container (width 150 mm), at the 
bottom of which a disk of radius R s — 45 mm rotates at an 
angular velocity Q ranging from 0.005 to 0.05 rpm, driven 
by a microstepping motor. 



Imaging of the suspension is performed using an index- 
matching technique, as introduced in 13.24-27. The 

PMMA 



as introduced in 13 24-2 
granular suspension consists of monodisperse 
spheres (Engineering Labs, diameter 4.6 mm ±0.1 mm 
placed in a mixture of Triton X-100, fluorescent dye (Nile 
Blue) and a few droplets of a 37% HC1 solution to tune 
the absorption spectrum of the dye. The fluid is prepared 
to closely match the refractive index of particles. The dif- 
ference in density between the fluid and the particles is 
A^llO ± 5 kg/m 3 . The fluids viscosity is 0.30 ± 0.05 
Pa.s. The split-bottom cell is filled with grains up to a 
fixed height H = 25± 2 mm; there is always a layer of fluid 
above the grains, so that the fluids surface tension does not 
play a role for the particle dynamics. The suspension is 
illuminated with a laser sheet (A=635 nm, power 30 mW) 
parallel to the bottom of the box (Fig. la) and mounted 
on a z-stage which allows the illumination of slices of the 
suspension at different heights (every 500 /im, starting 3 
mm from the bottom of the box). The thickness of the 
laser sheet is w 200 /xm. 

Image acquisition is done with a triggered 8 bit CCD 
camera operating at 10 Hz. Both the camera and laser 
sheet are mounted on vertical stages and can be trans- 
lated synchronously (velocity ~ 5 mm /sec) so that the 
probed slice of the suspension always stays in focus. Since 
the fluid is fluorescent but the particles are not, a single 
slice shows a collection of black disks (cross sections) on 
a white background (see Fig. lb). The contrast is suffi- 
cient to image the entire box, but best on the half closest 
to the laser, where we focus our analysis. We determine 
the particle trajectories using 3D tracking techniques de- 
veloped for confocal images of colloidal suspensions [28 



We follow the particles up to 1000 consecutive scans. We 
estimate the tracking error, obtained by tracking particles 
in a stationary packing, to be 1 /30 of the diameter in the 
plane of the sheet and 1/10 in the vertical direction. 

Flow profiles and local strain rates — From the 3D par- 
ticle tracking, we obtain a time sequence (maximum 1000 
steps) of particles positions fi (t) , where i denotes the la- 
bel of the approximately 5000 particle trajectories which 
can be reconstructed. Since the shear-rate is not homo- 
geneous in the split-bottom geometry, we divide the mea- 
surement area into bins. Given the azimuthal symmetry 
of the problem, these bins are chosen as squared tori (axis 
of revolution the vertical axis z) of size dji = dz=^ mm 
(rs 0.85 d). We determine the average velocity (v(r,z)) 
within each bin. The flow is orthoradial, as the trans- 
verse components v r and v z are negligible compared to 
vg — roj(r, z). 

The average angular velocity field oj(r, z) is independent 
of in the range used here 



of the ratio H/R s Wl 13 



13 



23 



As expected for this value 
the flow has a "trumpet 
like" shape, with shear bands arising from the edge of 
the rotating disc and propagating within the bulk of the 
suspension. We successfully fit cj(r,Zo)/Q for different 
Zq using error functions 12 13,23 (not shown). We use 
these fits to get an accurate estimate for the local strain 
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Fig. 2: (a) Examples of non dimensional mean squared displacements {As 2 ) /d 2 as function of time interval At. Q is fixed 
at 0.05 rpm, and we only plot data at z — 7 mm and z = 19 mm, for several values of r as indicated (see Fig. Id) — for 
these values, the local strain rate 7 varies over 2.5 decades. At late times, the curves all show a simple linear growth including 
diffusive behavior, with a diffusivity which clearly grows with the local strain rate, (b) When plotted as function of the local 
strain 7 := 7 At, mean squared displacements for the same value of z collapse, indicating that the dynamics are quasi-static. 
Families of collapsed curves are surprisingly distinct for different z. At late times, the curves all show a simple linear growth 
with strain, allowing the unambiguous definition of a dimensionless diffusion coefficient as D := (As 2 ) /d 2 /j. (c) By rescaling 
the local strain with the inverse of the dimensionless layer number z/d, the mean squared displacements for z — 7 mm and 
z = 19 mm can be scaled on top of each other, (d) Similarly, the values of D for all values of z and r (see Fig. Id) vary linearly 
with the inverse of the depth. The solid line is the prediction from Eq. 4. Inset: D as a function of -yd 2 , from which the pressure 
effect on diffusion is explicit. 



rate within each bin, assuming that the velocity gradient 
tensor has an unique non-zero component 



ire{r,Z ) = r/2 



du(r, gg) 
dr 



(1) 



The resulting local strain rate is plotted in Fig. Id. By 
varying the disk rotation speed from $7=0.005 rpm to 0.05 
rpm, we access local strain rates from 10~ 5 up to 10 -2 
a- 1 . 

Diffusivity — We have investigated the time evolution 
of the non-affme part of the trajectory : As^(At) = 
Ax^Ai) — (v^) At, where \i denotes the component r, 9 or 
z of cylindrical coordinates, and Ax^At) — x^(t + At) — 
x^(t). In the longitudinal (9) direction and for small values 



of r (see Fig. Id), the mean flow dominates, Taylor dis- 
persion plays an important role, and reliable data for the 
non-affine fluctuations are hard to obtain. We have thus 
focused on the z and r components of the fluctuations for 
large r, where the mean flow is slow (see Fig. Id). While 
all our data for the z-diffusion is consistent with the trends 
seen for the r-diffusion, the spatial resolution in z is less, 
so we focus our presentation on the radial diffusion. 

The normalized mean square displacements 
(As r (Af) 2 ) /g? 2 are plotted as function of At in Fig. 
2a, for two values of z, six values of r and a concomitant 
range of local strain rates which spans roughly two 
decades. At late times, the observed behavior is diffusive, 
and the mean squared displacement increases linearly 



P-3 



E. Wandersman et al. 



Q 
Q- 



1.00 


r_y 2 =43.6 a a): 






; _ y 2 =27.0 A 














10 




— Y 2 = 4.4 III ^ 




0.10 


r_y 2 =1.1 1 W 








c\j 






1 


0.01 










10" 1 




10' 



A ; A 2 1/2 

Ar/<Ar > 



10" 

Y 




Fig. 3: (a) PDF of radial displacements, for different values of the reduced kurtosis 72, showing the range of displacement 
distributions that we observe, (b) The radially averaged kurtosis decays with the 7.(0) The radially averaged kurtosis collapses 
when plotted as a function of the renormalized strain fd/z, and decays as a nontrivial powerlaw. 



with the time lag. We thus define a diffusion coefficient 
D as the limit for large At of the ratio (As r (At) 2 ) /(At). 
For the cases shown in Fig. 2a, this diffusion coefficient 
varies by two orders of magnitude. 

In the simple local, quasi-static picture, the mean square 
displacements would be determined by the local accumu- 
lated strain, 7 := jAt. Hence, one expects that these 
curves would collapse when plotted as function of the 
strain. Surprisingly, plotting (As,.(Ai) 2 ) /d 2 as function 
of 7 At does not achieve to collapse the data into a single 
curve, but rather we find two distinct curves, correspond- 
ing to data taken at the two different z positions (see Fig. 
2b). The diffusivity thus depends both on the strain rate 
and the vertical position in the sample. 

An alternative perspective on diffusivity starts from as- 
suming a Stokes- Einstein relation: D = kT / (Sirr/d) where 
D, rj and d are diffusion constant, effective viscosity and 
particle diameter, and kT refers to an effective tempera- 
ture 29 30 . For slow granular flows, shear stresses r are 
set by friction: r = /iP, where fi is the friction coefficient. 
As a consequence, the ratio of shear stress and strain rate 
(the effective viscosity), scales as 77 = nP/j, and substi- 
tuting this into the Stokes-Einstein relation, one obtains 
that 

D = (kT/3n^d)(i/P) , (2) 

so that the diffusion would be inverse proportional to the 
local granular pressure, which is set by the depth of the 
granular bed z (Eq. [2| . 

As a non-dimensional characteristic of the depth of the 
granular bed we take z/d — an integer that counts the 
number of grain layers to the free surface. To obtain 
data collapse consistent with the prediction D ~ P^ 1 , 
we can either multiply the mean squared displacements 
(As r (At) 2 ) with z/d, or multiply the local strain with 



d/z; as we explain below, we believe the latter captures 
the physics best. In Fig. 2c we show that all the mean 
squared displacements show good data collapse when plot- 
ted as function of jAtd/z. 

In Fig. 2d we combine our measured radial diffusivities 
for six values of r, five values of z and three global driving 
rates, and find that in good approximation: 



D ~ (jd 2 )(d/z) 



(3) 



with a proportionality constant of order one. If the diffu- 
sion coefficient is plotted as a function of the strain rate 
only, the diffusivities show a clear trend with z (see inset 
of Fig. 2d). 

Whereas there is mounting evidence for the nonlocal 
nature of slow granular flows [4 16 18; it appears that 
local plastic events affect the viscosity in regions further 
away. One might have expected that such nonlocal rhe- 
ology would go hand in hand with a nonlocal picture for 
diffusion, where the local diffusion constant is governed by 
the strain rate in a finite region around the point where 
the diffusion is probed. The good scaling collapse obtained 
with purely local metrics suggests that nonlocal effects in 
particle diffusion are weak or absent. 

Displacement Distributions — The rescaling of the dif- 
fusion constant with inverse pressure leaves open one ques- 
tion: is it fundamentally the strain scale or the displace- 
ment scale that is affected by the pressure? To answer 
this question, we have probed the probability distribution 
function (PDF) of the radial particle displacements over a 
range of lag times At, depths z/d and local strain rates 7. 

While for late lag times the distributions become Gaus- 
sian, for small At we see appreciable deviations, with very 
wide PDFs for the shortest lag times we can probe (See 
Fig. 3a). 
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We characterize the non-Gaussianity of the PDFs by 
the reduced kurtosis 72 = (Asf.) / (As 2 .) 2 — 3, and deter- 
mine the radial averages of the kurtosis, which thus only 
depends on z and 7. In Fig. 3b and 3c we compare plots 
of the kurtosis 72(2) for z — 6, 10, 14, 18 and 22 mm, as 
function of strain 7 and rescaled strain jd/z. Clearly, the 
rescaled strain performs a better rescaling of the kurtosis, 
even though the curve for z — 22, i.e. close to the bot- 
tom, appears to deviate from the other curves. Together 
with our rescaling of D, this provides strong evidence for 
a pressure dependent characteristic strain jAtd/z. 

We finally note that on short time scales, the parti- 
cle displacement distributions appear to have fat tails. 
We note that the kurtosis shows clean power law scal- 
ing with the rescaled strain, with an exponent of approx- 
imately -0.7. Such power law decay of the kurtosis may 
happen when the underlying stochastic process is heavy- 
tailed [32], and is a signature that the distribution of in- 
dividual jumps is of power law form. 

Discussion — Our results show that particle displace- 
ments in sheared granular flow are diffusive on long time 
scales, with a diffusivity proportional to the shear rate 
[l][9] and a surprisingly inverse proportionality to depth. 

A natural question concerns the microscopic origin of 
the diffusion scale in slow granular flows. In a Langevin 
equation type of approach, the order of magnitude of the 
diffusion coefficient is D ~ where F is the stochastic 
force acting on the particles and 77 = the effective drag 
coefficient. 

Considering the forces (gravitation, viscous, friction, in- 
ertial) and the velocity scales (say U < R S Q) in the present 
study, we see that the viscous forces F v ~ ZirrifdU (where 
rjf is the fluid viscosity) and inertial forces Fi ~ ^ir - ^ - 
are negligible with respect to gravitational F g — ^^-Apg 
and frictional forces Ff ~ pPd 2 . In other words, the 
Reynolds number is very small Re = Fi/F v ~ 10 -4 and 
the Stokes number is large St — F g /F v ~ 10. We are in 
the "free fall" regime described by Courrech du Pont et. 
al. in 
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indeed, we expect our data to be relevant for 
slow dry flows also. 

This leaves frictional and gravitational forces as scales 
that set the diffusion. Clearly, friction sets the drag, but 
it is not obvious whether friction or gravity provides the 
dominant stochastic force. If friction would dominate, 
both the driving and damping would be set by it, and 
we would expect that the diffusion coefficient would be 
independent of depth. In contrast, the combination of a 
gravitational drive and frictional damping gives the ob- 
served depth dependence of the diffusion as we show now. 
Assuming that gravity sets the driving force for diffusion, 
we find that: 
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D = -jd 

6 pp 



jd 2 - 
z 



(4) 



since p sa Apgz and p s» 0.5 13 . This expression is plot 



ted as a solid line in Fig 2d. The agreement is reasonable, 



supporting the idea of gravity as a driving force for gran- 
ular diffusion in slow flows. 
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